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Abstract. The method of elevation classes, in which the ice 
surface model is run at multiple elevations within each grid 
cell, has proven to be a useful way for a low-resolution atmo- 
sphere inside a general circulation model (GCM) to produce 
high-resolution downscaled surface mass balance fields for 
use in one-way studies coupling atmospheres and ice fiow 
models. Past uses of elevation classes have failed to conserve 
mass and energy because the transformation used to regrid 
to the atmosphere was inconsistent with the transformation 
used to downscale to the ice model. This would cause prob- 
lems for two-way coupling. 

A strategy that resolves this conservation issue has been 
designed and is presented here. The approach identifies three 
grids between which data must be regridded and five trans- 
formations between those grids required by a typical cou- 
pled atmosphere-ice fiow model. This paper develops a the- 
oretical framework for the problem and shows how each of 
these transformations may be achieved in a consistent, con- 
servative manner. These transformations are implemented in 
Glint2, a library used to couple atmosphere models with ice 
models. Source code and documentation are available for 
download. Confounding real-world issues are discussed, in- 
cluding the use of projections for ice modeling, how to han- 
dle dynamically changing ice geometry, and modifications 
required for finite element ice models. 


1 Introduction 

Many questions still surround the issues of how ice sheets 
respond to climate forcing and how those changes will af- 
fect sea level, regional and global climate. Recent observa- 
tions have shown accelerating mass loss from the Greenland 
and Antarctic ice sheets (Vaughan et al., 2013), adding ur- 
gency to these questions. Although general circulation mod- 
els (GCMs) are able to project changes in ice sheet surface 
mass balance, projection of changes in ice sheet mass due to 
ice dynamics with GCMs remains a challenge (Church et al., 
2013). A number of climate modeling groups are addressing 
these deficiencies by adding dynamic ice fiow effects to ex- 
isting GCMs: GISS ModelE (Schmidt et al., 2006), CESM 
(Hurrell et al., 2013), HadGEM2 (Collins et al., 2011), etc. 
This is done by coupling the GCM with an existing ice flow 
model, such as Glimmer-CISM (Rutt et al., 2009), BISI- 
CLES (Comford et al., 2013), ISSM (Larour et al., 2012), 
PISM (Bueler and Brown, 2009), etc. A full understanding 
of the long-term evolution of an ice sheet within a coupled 
climate system requires coupling with the ocean as well as at- 
mosphere. Surface runoff, ocean cavity circulation and salin- 
ity gradient effects are all important. In this paper, we focus 
only on coupling with the atmosphere. 

One can distinguish between one-way and two-way cou- 
pling. In one-way coupling, the GCM is used to develop sur- 
face mass balance (SMB) and temperature fields, which are 
then used to drive the ice fiow model off-line. This process 
misses effects caused by feedbacks from the ice sheet to the 
rest of the Earth system: for example, decreased ice sheet 
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albedo (Qu and Hall, 2006) or lowered atmosphere orogra- 
phy (Ridley et al., 2005). Past studies with one-way coupling 
have yielded useful insight into the future of present-day ice 
sheets; examples include Huybrechts (1994), Greve (2000), 
Stone et al. (2010), Bindschadler et al. (2013), Lipscomb 
et al. (2013), Nowicki et al. (2013a), Nowicki et al. (2013b) 
and Goelzer et al. (2013). However, ice sheet feedbacks are 
expected to be increasingly important for simulations of the 
long-term evolution of ice sheets and the climate associated 
with them. This kind of feedback probably plays a significant 
role in many events in the paleorecord (Dansgaard-Oeschger 
events, Heinrich events, the Younger Dry as). 

Two-way coupling strategies address these issues by al- 
lowing the atmosphere to be influenced by changes in the 
ice sheet elevation, extent and albedo over time. We distin- 
guish between loose and tight two-way coupling. Loose two- 
way coupling involves running a series of GCM simulations 
with different ice sheet configurations, each based on the re- 
sult of the previous. Each GCM run is a separate simulation, 
without continuity of mass or energy between runs. Stud- 
ies with loose two-way coupling have yielded insight into 
future equilibrium states for ice sheets and climate (Ridley 
et al., 2005, 2010). However, ice sheets change configura- 
tion in these runs without accompanying mass and energy 
fluxes required to make those changes happen. This is equiv- 
alent to applying an unknown impulse forcing to the ice sheet 
each time it is changed. Although the coupled ice sheet might 
eventually reach the correct equilibrium state, the transients 
involved in achieving that equilibrium will be suspect. Un- 
fortunately, results relevant to human society all require an 
understanding of the transients. For successful simulation of 
transients, we turn to tight two-way coupling. It involves run- 
ning the ice flow model, step by step, along with the rest of 
the GCM - while conserving mass and energy along the way. 
Attention to conservation is required, since the GCM is sim- 
ulating a more nearly closed system that could run for a long 
time. 

When one couples dynamic ice flow models with GCM 
atmospheres, two models operating on different grids and 
timescales must communicate: ice flow models operate at 
low frequency on a high-resolution grid with local projec- 
tion, while GCM atmospheres operate at high frequency on 
low-resolution global grids. A number of issues arise due 
to this mismatch, including how one creates high-resolution 
surface mass balance fields from low-resolution GCM in- 
put. Elevation classes address the latter issue. They were first 
introduced for precipitation downscaling in a GCM by Le- 
ung and Ghan (1998) and later applied to one-way coupling 
from GCM atmosphere to ice flow models by Lipscomb et al. 
(2013). The key insight is that mass and energy fluxes be- 
tween the atmosphere and an ice sheet vary approximately 
by elevation within a local region. 

When using elevation classes, a third grid is introduced, 
the elevation grid. This allows the GCM to compute surface 
fields at a variety of elevations within each atmosphere grid 


cell, not just the elevation seen by the atmosphere. A high- 
resolution surface mass balance is produced on the ice grid 
by first computing SMB on the elevation grid, and then us- 
ing a vertical interpolation scheme to produce SMB on the 
ice grid. This method of interpolation produces surprisingly 
good results: although it cannot capture certain localized ef- 
fects (e.g., wind direction), it has been shown to allow GCMs 
to produce surface mass balance fields approaching the qual- 
ity of those produced by regional climate models (Vizcaino 
et al.,2013). 

Inside the GCM, SMBs computed on the elevation grid 
must be regridded to the atmosphere grid as well as the ice 
grid. In order to maintain conservation in a tight two-way 
coupled system, it is essential that the set of regridding op- 
eration chosen is self-consistent: that is, if a flux field on the 
elevation grid is regridded simultaneously to the atmosphere 
and ice grid, then the total amount of flux represented by 
the resulting two fields should be the same. For conservation 
purposes, the specifics of these two transformations are not 
important, as long as they are consistent with each other. Past 
efforts at one-way coupling have defined these two transfor- 
mations in ways that each make intuitive sense, but are not 
consistent with each other. This is not a problem for one-way 
coupling, but it would cause conservation problems for two- 
way coupling. 

This paper develops the concept of the elevation grid on 
which the ice surface model runs, and then derives a set of 
conservative regridding transformations between the atmo- 
sphere, elevation and ice grids. The coupled processes un- 
der consideration along with the transformations required for 
tight two-way coupling are introduced in Sect. 2. Section 3 
focuses on the elevation grid, while the grid fundamentals 
necessary for the two-way coupling are presented in Sect. 4. 
Section 5 deals with the use of projection and associated is- 
sues encountered when bridging between the spherical ge- 
ometry of GCMs and Cartesian ice sheet models. We show 
how to choose and implement the transformations in Sects. 6 
through 9 and work through realistic examples of these trans- 
formations in Sects. 10 through 12. We touch on a number of 
extra “wrinkles” in the real-world problem: procedures to use 
when the elevation grid is based on a horizontal grid other 
than the atmosphere grid in Sect. 13, and regridding proce- 
dures required when ice elevations, ice extent or elevation 
grid change in the simulation in Sect. 14. Finally, we present 
in Sect. 15a library, Glint2, that can be used to tightly couple 
GCMs and ice sheet models. 


2 Coupled processes 

The coupled atmosphere-ice system involves three models 
interacting with each other: an atmosphere model, an ice 
flow model and an ice surface model situated between them 
(Fig. 1). The atmosphere model is coupled with the ice sur- 
face model, which tracks the top few meters of ice. Processes 
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Fig. 1. Configuration of the three models, and the separate spa- 
tial domains they oeeupy. The atmosphere model (A) is eoupled 
to the iee surfaee model (F), whieh traeks the top few meters of 
iee. The iee surfaee model (F) is then eoupled to the dynamie iee 
flow model (I). The dynamie iee flow model, operating at long time 
steps, is insulated from high-frequeney surfaee proeesses in the iee- 
atmosphere interaetion. 


modeled here include a full surface mass-energy balance 
computation, snow-fim compaction, albedo effects, meltwa- 
ter percolation, runoff, refreezing, etc. The bottom of the ice 
surface model is coupled to the ice flow model. 

2.1 Ice surface model 

In order to couple an ice sheet to a GCM atmosphere, SMB 
on the ice sheet must be calculated from GCM outputs. 
Some one-way coupled studies have used temperature index 
schemes (Huybrechts, 1994; Greve, 2000; Stone et al., 2010; 
Bindschadler et al., 2013; Nowicki et al., 2013a, b; Goelzer 
et al., 2013): the mean surface temperature over each cou- 
pling time step (typically one month or year) is used to com- 
pute SMB, and both are passed to the ice flow model. It is 
hard to see how energy can be conserved with such a scheme. 
The problem is that the atmosphere must run for many time 
steps before atmosphere-ice sheet energy fluxes are com- 
puted on a coupling time step. The computed flux on the ice 
sheet will not be the same as that sum of fluxes seen by the 
atmosphere over the previous coupling time step - and it is 
“too late” to go back and change the atmosphere to match. 

For this reason, a full energy balance scheme, computed 
each atmosphere time step (typically one hour) and inte- 
grated over the coupling time step, is considered essential 
in a GCM setting. Energy flux between atmosphere and ice 
sheet follows diurnal and seasonal cycles, making positive 
and negative contributions to the integrated flux. It is impor- 
tant that energy flux is computed at a small enough time step 


to capture these effects accurately; typically, the atmosphere 
time step is suflflcient. 

While ice sheets move slowly and can be modeled with 
long time steps, the modeling of the surface of the ice sheet 
requires short time steps. This is accomplished by introduc- 
ing the ice surface model as a third model, sitting in between 
the ice flow model and atmosphere model. The top of the ice 
surface model couples with the atmosphere model every at- 
mosphere time step, whereas the bottom couples with the ice 
sheet model every coupling time step. The ice surface model 
needs to be thick enough so there is little variation in tem- 
perature at the bottom: 15 m is sufiicient. This ensures that 
the heat flux with the ice flow model will be small and low 
frequency. 

The ice surface model can serve an additional purpose of 
modeling the great variety of surface processes that may be 
relevant to the long-term evolution of ice sheets: snow-flm 
compaction, surface runoff/drainage networks, water per- 
colation, refreeze, albedo effects and wind-blown snow, to 
name a few. 

The use of an ice surface model solves some important 
problems, but it also introduces nonphysical elements into 
the model. Ice contained in the ice surface model does not 
advect along with the ice flow model, nor does it contribute 
to the stress fleld of the ice below it in the ice flow model. 
In both cases, we expect the relative error to be small: the 
top 15 m of snow/flm contains less than I % of the mass of 
a 1000 m thick ice sheet. There may be ways to flx these 
problems - however, there is no need to make the ice surface 
model thicker than it needs to be. Numerous studies have 
shown that ice sheets below about 1 5 m are fully insulated 
from surface weather and seasonal cycles (Zagorodnov et al., 
2006). 

2.2 Three models, three grids 

Each of the three models in the coupled system runs on its 
own grid. The atmosphere is run on the atmosphere grid (A) 
and the ice flow model is run on the ice grid (/). The ice sur- 
face model is run on the elevation grid (E), which is based 
on the atmosphere grid (see Sect. 3). All three grids are two- 
dimensional, in the sense that they are used to construct two- 
dimensional functions f(x,y) over the domain. Regridding 
operations are needed to pass mass and energy fluxes be- 
tween the models. 

It is important to keep in mind the relative size of the 
three grids. We set up a test using the GISS 2° x 2^° atmo- 
sphere grid (Schmidt et al., 2006), overlapping with the 5 km 
grid from SeaRISE (Sea-level Response to Ice Sheet Evolu- 
tion; Bindschadler et al., 2013). We used 40 elevation points, 
spaced every 100 m from Om to 4000 m. In this case, A had 
146 grid cells, I had 66 906 and E had 1829. These numbers 
only account for grid cells involved with the Greenland ice 
sheet. In general, the ice grid will be finest, the atmosphere 
grid coarsest and the elevation grid in the middle. 
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Fig. 2. Data flow (blue arrows) for the eoupling between atmo- 
sphere, iee surfaee and dynamie iee flow models (boxes). Sinee the 
three models run on different grids, regridding operations (ovals) 
are required at eaeh step. Figure 21 shows the inputs required to 
eompute these regridding operations. 

Time frequency mismatches are another issue. The atmo- 
sphere runs at high frequency, each time step being typically 
1 h. On the other hand, the atmosphere and ice flow models 
are coupled at much lower frequency, typically one month 
or any other time period. We call these two time steps the 
atmosphere time step and ice coupling time step. 

2.3 Fully coupled system 

Figure 2 shows the data flow of the fully coupled system. 
Steps of the data flow are organized based on their frequency: 
the top circle of steps runs at the same frequency as the GCM 
atmosphere, while the bottom circle runs at the ice coupling 
frequency - typically one month or more. We describe the 
steps involved in coupling the GISS ModelE (Schmidt et al., 
2006) with PISM (Bueler and Brown, 2009); however, these 
steps are general for any GCM or ice flow model. We now 
trace through the data flow on a typical coupled run. 

2.4 Atmosphere time step 

When the atmosphere runs, it produces a set of flelds on 
the atmosphere grid that affect the processes in the ice sur- 
face model: for example, downwelling long-wave radiation. 


A Atmosphere grid, projected to Cartesian plane 
A' Original atmosphere grid on the sphere 
I Ice grid 

G Interpolation grid (e.g., ice or exchange grid) 

E Elevation grid 
M Interpolation matrix: E ^ G 
R Area- weighted remapping matrix: G ^ A 

A “Repeat” transformation matrix: A^ E 
X Area- weighted remapping matrix: G ^ I 

X' Area- weighted remapping matrix: I ^ G 

P Diagonal scaling matrix: A! ^ A 
LO Function defined to be constant within grid cells 
LI Function defined by piecewise linear 
interpolation between grid points 

Fig. 3. Definition of symbols used throughout the paper. See Ap- 
pendix A for notational conventions. 


downwelling shortwave radiation, precipitation, etc. These 
fields must be regridded from A to E. We denote the set of 
fields to be regridded with a capital letter, ; the result of 
the regridding is R^ (see Appendix A and Fig. 3 for nota- 
tional conventions). 

The ice surface model is run at the same frequency as the 
atmosphere. Among its outputs are mass and energy fluxes 
with the atmosphere: evaporation, sublimation, upwelling 
long-wave radiation, latent heat release, etc. These are repre- 
sented by — F^ and are regridded to — F"^ before being passed 
back to the atmosphere on the next time step. 

The ice surface model also produces fluxes in the ice flow 
model; these are described immediately below. 

2.5 Coupling time step 

On each atmosphere time step, relevant flux outputs from the 

- E 

ice surface model are accumulated as F for future coupling 

- £ 

with the ice flow model. They are named F because these 
fluxes are in general equal and opposite to fluxes sent to the 
atmosphere. Every coupling time step - about once a month 

- E -I 

- the accumulated F is regridded to the ice grid (F ) and 
passed to the ice flow model. 

The ice flow model produces changes in ice surface topog- 
raphy and extent, as well as a small energy flux between the 
ice flow and ice surface models (together, we call these D^). 
Changes in ice topography and extent are regridded to the 
atmosphere grid (D^) and used to adjust the atmosphere’s 
orography. The energy flux is regridded to the elevation grid 
(D^) and applied to the ice surface model. 

2.6 Regridding requirements 

We see that the fully coupled system requires the use of five 
regridding operations at various points in its run (Fig. 2). 
These regridding operations must be conservative, in the 
sense that none of them can change the integral of the field 
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on the domain as a whole. In fact, we would like to impose a 
stronger conservation condition so that values are conserved 
within each atmosphere grid cell. 

We develop these five regridding operations in the sec- 
tions below. We begin by discussing the method of elevation 
classes in a general manner (Sect. 3) and then move on to 
developing basics of conservative regridding (Sect. 4). Be- 
cause conservation is defined in terms of integration over ar- 
eas, part of our discussion involves a definition of how to 
integrate functions on the elevation grid E. Finally, we show 
how the required set of operations can be constructed in a 
fully self-consistent manner. 

Throughout, we assume that the atmosphere grid A has 
LO parameterization: functions on A are represented by their 
mean value within each grid cell. This becomes our basis for 
conservation: all regridding operations are conservative over 
every atmosphere grid cell. We do not require any specific pa- 
rameterization for the ice grid, but we consider cases for LO 
(piecewise constant) and LI (piecewise linear) parameteriza- 
tions, the latter which are commonly used in finite element 
ice models. 

3 Elevation points 

The method of elevation classes, when applied to the 
atmosphere-ice coupling problem, involves running the ice 
surface model at one or more elevations in each atmosphere 
grid cell (Lipscomb et al., 2013). Temperature, pressure and 
precipitation are extrapolated to a set of elevations within the 
grid cell. They are then used in an ice surface model to com- 
pute a full surface mass and energy balance at each elevation. 
The result is a set of “what-if ’ scenarios, giving an estimate 
of what the fiuxes would have been, had the ice surface of the 
grid cell been at the given elevation - rather than the eleva- 
tion seen by the atmosphere model. 

The modeler must choose which elevations to use for each 
atmosphere grid cell. The simplest approach is to use a fixed 
set of elevations across all grid cells - for example, every 
100 m from 0 m to 4000 m. 

However elevations are chosen, the result is a new “grid” 
- the elevation grid - on which the ice surface model is run 
and surface fiuxes are generated. The elevation grid is derived 
from the atmosphere grid, in the sense that each elevation 
grid “cell” (or elevation point) is associated with one parent 
atmosphere grid cell. If elevation point Ej is associated with 
atmosphere grid cell A/, we write Ej g At. 

Once the GCM has computed a conserved quantity on the 
elevation grid, those values can be used to develop a relation 
between elevation and SMB within each atmosphere grid 
cell. Suppose we have computed a fiux field f^ \ SMB, for 
example. For an atmosphere grid cell A/, consider the com- 
ponents of the fiux field that are related to the enclosing 
atmosphere grid cell A/. Or more formally, consider for 
all j such that Ej g A/. We can think of these component 



Fig. 4. Interpolated surfaee mass balanee (SMB) funetion within 
one atmosphere grid eell. These values eome from one month (July) 
of a run of GISS ModelE with elevation points. The vertieal dashed 
line indieates the mean elevation of the grid eell, “seen” by the at- 
mosphere. Dots and squares represent the results of extrapolated 
SMB eomputations at other elevations: dots represent elevations 
found within the grid eell, whereas squares represent elevations out- 
side the eelFs range. SMB values at intermediate elevations are in- 
terpolated. 

values as samples of a 1 -D function relating elevation to fiux 
(SMB), within the localized region of At. By interpolating 
between those points using standard methods, one can con- 
struct a continuous function relating fiux to elevation within 
Ai . As long as enough points are used and the function being 
interpolated is smooth enough, this procedure will yield an 
arbitrarily accurate representation of the “true” function. 

In fact, the functions we typically wish to interpolate - sur- 
face mass and energy balance averaged over about a month - 
are quite smooth as functions of elevation (Fig. 4). In the face 
of spatially invariant precipitation, one would expect SMB 
to be constant to first order above the equilibrium line alti- 
tude (ELA), and to decrease linearly with elevation below 
the ELA. Near the ELA, one would expect a smooth tran- 
sition because the ELA goes up and down over a month of 
diurnal cycles. This is in agreement with experimental work, 
which has shown SMB below the ELA to vary linearly with 
temperature (Braithwaite, 1981; Box et al., 2013). 

How many elevation points are required to properly re- 
solve the elevation-fiux relationship? This depends on the in- 
terpolation scheme: higher order schemes will require fewer 
points than piecewise linear interpolation. But the general 
shape of the function - two straight lines connected by a 
curve near the ELA - implies that not many points are 
needed, except for near the ELA. 

If one knows where the ELA is on every atmosphere grid 
cell, then this is a useful guide in selecting elevation points. 
But if one is studying ice sheets in a changing climate, then 
the ELA will be expected to move over time. It is possible 
to adaptively move elevation points as the ELA moves. But 
it is simpler just to use many points everywhere. In our tests, 
we have used 40 points at 100 m spacing, which is proba- 
bly more than sufihcient for coupled ice sheet simulations. 
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We have not done a careful study of the optimal number of 
elevation points. 

The method of elevation classes - or elevation points - 
provides a way to construct an interpolated relation between 
elevation and surface mass-energy flux within each atmo- 
sphere grid cell. Additional choices need to be made in or- 
der to produce fully downscaled flux fields on the ice grid 
(Sect. 6). We are now almost ready to define the transforma- 
tions posited in Sect. 3. But first, we pause for some ground- 
work on grid fundamentals in Sect. 4 and projection issues in 
Sect. 5. 

4 Grid fundamentals 


as needed. How one integrates the basis functions depends 
on the nature of those basis functions. We give specifics for 
LO and LI grids in Appendices B and C, respectively. 

4.2 Comparison across grids 

In regridding applications, we need to compare fields across 
different grids. They cannot be tested for simple point-by- 
point equality, due to differences in grid structure. Instead, 
we compare by integrating two fields over the same region or 
set of regions. 

If we have two fields and on two different grids, 
we say that is equivalent to on region B, ox 

if 


Numerical models represent continuous fields with linear 
combinations of a finite set G of basis functions - which we 
call a parameterization (Appendix A describes our mathe- 
matical notation). For example, suppose that G uses the basis 
functions g(x, y) = [gi{x, y), . . . , gn(x, y)], where g(x,y) 
is the vector of all the basis functions. Suppose we have an 
^-dimensional vector with components /.^. That vector 
represents the function f^{x, y) formed by a weighted sum 
of the basis functions: 

P(x,y) = P ■g(x,y). (1) 

A wide variety of basis function sets are used for different 
problems. Most commonly in climate modeling, repre- 
sents the mean value of f^(x,y) within some well-defined 
bounded region, which we call a grid cell. From a conserva- 
tion point of view, a function f^(x,y) with LO parameteri- 
zation may be taken to be constant within each grid cell, with 
discrete “jumps” from one grid cell to the next. In this case, 
the basis function gt (x , y) for grid cell i is equal to 1 inside 
the cell and 0 outside. 

LO parameterizations are widely used in climate and fi- 
nite difference ice flow models. Finite element ice flow mod- 
els might use LI or even higher order parameterizations 
(Zienkiewicz et al., 2013), and they use the term “mesh” to 
describe the geometry of their basis functions. In this paper, 
we use the term “grid” to refer to both the vector space in 
which lives and its associated parameterization. 

4.1 Integration on grids 

Because we are working with conserved quantities, it is es- 
sential that we can integrate functions over any well-defined 
region B. Because a parameterization defines f^(x,y) on 
every point, this integral is well-defined. Substituting from 
Eq. (1) and using linearity, we get 

j P(x,y)dA = p-jg(x,y)dA. (2) 

B B 

If we wish to integrate over B repeatedly, we can pre- 
compute J^g(x,y)dA, and then take dot products with 


P ■ j g(x, 3')dA = ■ j h(x, 3?)dA. (3) 

B B 

Suppose we wish to compare and over an entire 
domain? If we have a set of regions A = {A\, . . . , An} tiling 
that domain, then we can say the two fields are equivalent 
on A if they are equivalent on all of Ai , . . . , A„. Note that A 
could be an LO grid, or simply a set of regions on the domain. 
If we have a regridding transformation such that f^ = 
*^(/^)? then we say that is conservative on A if 


4.3 Area-weighted remapping 


Suppose we have two grids G and H, with H being LO pa- 
rameterized - for example, G might be an ice grid and H 
an atmosphere grid. Suppose we have a held /^, which we 
wish to regrid in a conservative manner to result in the held 
f^. One common way to do this is to compute each com- 
ponent f-^ based on the value of f^(x,y) integrated over 
the area covered by the /th grid cell in the grid H. Or, more 
formally. 




1 

m 


\ 

P ■ f g{x,y)dA 

Hi / 


( 4 ) 


By construction, this transformation is conservative on H. 
Note that and will generally not be equivalent over 
other regions other than grid cells m H - for example, over 
grid cells in G (if G happens to be LO parameterized). 

Area-weighted remapping is closely related to our defi- 
nition of equivalence above. If we have two fields on two 
different grids, equivalence on H can be tested by regrid- 
ding both fields to H and comparing the resulting vectors. 
The transformation from G to H is linear and thus repre- 
sentable by a matrix. It is variously called area-weighted 
remapping or conservative regridding (Ramshaw, 1985). The 
matrix is computed by integrating the source basis functions 
g\{x,y), . . . , gn{x,y) over every grid cell in the destination 
grid. 
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Fig. 5. An exchange grid, obtained by overlapping a sample atmo- 
sphere grid and ice grid (not to scale). Each resulting grid cell, ir- 
regularly shaped, overlaps with exactly one atmosphere grid cell 
and exactly one ice grid cell. 


4.4 Interpolation grid 

We are now almost ready to discuss our implementations of 
the transformations E ^ 1 and E ^ A (see Appendix A and 
Fig. 3 for mathematical notation). But first, we must intro- 
duce the interpolation grid G, which is used to rigorously 
define basis functions on the elevation grid E. The user has a 
choice of how G is to be chosen. 

In general, G is chosen simply as the ice grid /. If / is 
LO-parameterized, we might instead choose G to be the ex- 
change grid between A and I (Fig. 5): the LO grid whose 
grid cell outlines are formed by the intersection of grid cells 
in A and I (Balaji et al., 2006). The exchange grid is a use- 
ful choice for G because every exchange grid cell overlaps at 
most one atmosphere grid cell. This choice has its pros and 
cons, which we explore in Sect. 9. Either way, the interpola- 
tion grid G is similar to the ice grid I and can be thought of 
as an ice grid proxy in most cases. 


5 Projection issues 

Whether the source grid is LO or LI, area-weighted remap- 
ping algorithms need to find the intersection of grid cell out- 
lines from two grids. Technically, this is only possible if the 
two grids exist on the same surface. In our case, atmosphere 
models exist on the surface of a sphere, whereas ice flow 
models work on a Cartesian plane. Unless an ice flow model 
formulated in spherical coordinates is used, the exchange 
grid between an atmosphere grid and an ice grid cannot be 
directly computed. 

This problem is solved using a map projection (Snyder, 
1987). We let A' be the original atmosphere grid and then let 
A be the projected atmosphere grid - projected to the Carte- 
sian plane using the chosen projection. Regridding compu- 
tations described in this paper are made to/from A, the pro- 
jected grid. 

In general, the area of a grid cell can change when it 
is projected. Changes in area due to the use of a projec- 
tion that does not preserve areas are called projection error 


(Lauritzen andNair, 2008). Even if an equal area projection 
is used, in practice grid cells still experience small changes 
in area. This is because projected grid cells have complex 
shapes, but the algorithms in Appendix D are only able to ap- 
proximate these shapes with polygons. Changes in area due 
to polygonal approximation of curves, or numerical artifacts 
in the methods used, are called geometric error (Lauritzen 
and Nair, 2008). 

Whether a change in grid cell area arises from projection 
or geometric error, data must be rescaled to maintain conser- 
vation when transforming a Add to f^ \ 


fi^ = 



( 5 ) 


Note that this transformation is diagonal and invertible. 
Because this is simply a rescaling, the regridding schemes 
presented in this paper are not operationally affected by pro- 
jection issues. 

Although this rescaling scheme may be used to address 
any change in grid cell area, it is preferable to eliminate or 
minimize these changes to begin with. To this end, we con- 
sider projection and geometric error separately. 


5.1 Projection error 


The stereographic projection used in SeaRISE (Bindschadler 
et al., 2013) can change the area of grid cells by up to 10 %. 
If a projected grid cell is 10 % smaller then the original, then 
that means that 1 m of accumulation in the GCM will turn 
into 1 . 1 m in the ice flow model. This could cause significant 
discrepancies in dynamic ice flow. Because the projection in 
SeaRISE is constructed to minimize errors in a band along 
the TUN parallel, we would not expect the magnitude of 
projection error to change significantly if an oblique stereo- 
graphic projection were used instead. 

One solution to the problem of projection error would be 
to use an equal area projection - the Lambert equal area 
projection, for example. However, equal area projections do 
not preserve angles, which distorts the ice dynamics at a lo- 
cal scale. We conclude that nonphysical distortions happen 
whether an area-preserving or angle-preserving projection is 
used. It is not yet clear which choice gives better results in 
the end. 

One innovative approach to this problem is to use an angle- 
preserving projection, but to allow ice flow model grid cells 
to vary in size, based on the local distortion of space intro- 
duced by the projection. Parameters dx and dy are kept for 
each grid cell. These grid parameters are included in the ice 
flow model equations, thereby accounting for shape and area 
distortion caused by the projection in a physically meaning- 
ful way (Pollard and DeConto, 2012). 

5.2 Geometric error 


Geometric errors are typically three orders of magnitude 
smaller than projection errors. In theory, they may be 
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Fig. 6. An example of geometrie error, in whieh grid eells ehange 
size due to polygonal approximation. This map, made using a Lam- 
bert equal area projeetion eentered on the North Pole, shows a set of 
latitude-longitude grid eells - a kind eommonly used in atmosphere 
models. Solid blue lines show polygonal approximations, while dot- 
ted red shows the aetual grid eells on the sphere. Note that all grid 
eells (in this loeal map) shrink when approximated. 


arbitrarily reduced further by increasing the number of sides 
used in the approximating polygons - the parameter n spec- 
ifies the number of line segments used to approximate each 
side of the original grid cell in its projected state. But this 
method has its limit in practical terms. 

Consider a typical latitude-longitude grid on the sphere 
with a Lambert equal area projection centered at the North 
Pole (Fig. 6). In this case, the area of a polygonal approxi- 
mation will always be smaller than the area of the actual grid 
cell: we end up inscribing polygons inside of circles. The 
geometric error will depend on the number of sides of the in- 
scribed polygon, which in this case depends on the number 
of grid cells in the circle, and on n. 

This is the method used by Archimedes to approximate 
the value of tt in ~ 250 BCE (Heath and Archimedes, 1897, 
p. 91). Unfortunately, it converges only quadratically, as 
0{\/n^). Meanwhile, memory use to store all those poly- 
gons goes up by 0(n). Memory and time requirements will 
therefore be exponential in terms of the number of digits of 
accuracy required: each additional digit of accuracy will re- 
quire an increase in the number of sides by a factor of ^lO. It 
is therefore not practical to make geometric error arbitrarily 
small by increasing n. In our experience, a value oi n = 2 
yields geometric error of approximately 10“^. We believe 
n = 2 offers an acceptable trade-off between efficiency and 
accuracy. 


6 Transformation: elevation to ice/interpolation grid 

We can now derive the first of the five transformations 
posited in Fig. 2. We begin with the transformation E ^ G 
from the elevation grid to interpolation grid, because this is 
the transformation that ice modelers are most concerned with 
when seeking high-quality downscaled flux fields. In various 
contexts, this transformation might be referred to as an inter- 
polation, downscaling or regridding operation. 

In Sect. 3, we showed how values at elevation points 
may be used to construct a flux-elevation relationship fi (z) 
in each atmosphere grid cell A/. These relationships may 
then be used to construct on the interpolation grid. The 
resulting transformation for E ^ G potentially involves hor- 
izontal and vertical interpolation, of which the modeler has 
considerable choice. 

We present here three methods of interpolation from 
E to G: z interpolation (Sect. 6.1), bilinear interpolation 
(Sect. 6.2) and elevation class interpolation (Sect. 6.3). All 
three methods assume that G is LO parameterized and can be 
used for the basis of a conservative coupling. And since they 
are all linear, they can all be represented by a (sparse) matrix, 
which we will call M. The methods are extended to the LI 
case in Sect. 6.4 

6.1 Z interpolation 

Z interpolation constructs a flux field on the interpola- 
tion grid through direct application of the flux-elevation re- 
lationship fi (z) derived in Sect. 3 for atmosphere grid cell 
Ai . Supposing G is LO parameterized, consider an interpola- 
tion grid cell Gj, wholly contained in one atmosphere grid 
cell Ai, with mean elevation Zj. We would set the value of 
that grid cell to = f (zj). If Gj intersects more than one 
atmosphere grid cell, the same procedure is followed for each 
atmosphere grid cell Ai it intersects - and the results are 
summed together, area-weighted by | Gy H A/ 1 . 

We call this z interpolation. As long as standard interpo- 
lation methods are used to construct //(z), z interpolation 
defines a linear function from to and can be rep- 
resented by a matrix. An example of the result is shown in 
Fig. 7a. In general, z interpolation produces SMB fields that 
vary smoothly within each atmosphere grid cell but that con- 
tain discontinuities between grid cells. 

6.2 Bilinear interpolation 

There is some concern that discontinuities created by z in- 
terpolation could cause problems as an input to an ice flow 
model. In that case, a bilinear interpolation step may be used 
to create a smooth field, as in Lipscomb et al. (2013). This 
scheme sets the value of interpolation grid cell Gy equal to 
a linear combination of //(z), fk{z), fi(z) and fm(z), where 
atmosphere grid cells A/, A^, A/ and A^ are the four grid 
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- 20,00 0.00 10.00 - 20:00 0.00 10-00 
SMB (kg/day*m^2) SMB (kg/day*m^2) 

(a) (b) 


Fig. 7. Results of (a) z interpolation and (b) bilinear interpolation, 
to develop a downsealed field on the iee grid from an elevation point 
field. Bilinear interpolation eliminates the diseontinuities present at 
atmosphere grid eell boundaries. However, total SMB is ehanged, 
partieularly for loealized events sueh as snowstorms, for example in 
the northwest. This eould introduee biases in the long-term evolu- 
tion of the iee sheet, even though both sehemes ean be made to eon- 
serve mass and energy. This figure is for demonstration purposes 
only. See Seets. 10 and 1 1 for thorough regridding examples. 

cells with centers closest to the center of Gj (Fig. 8). Results 
are shown in Fig. 7b. 

Bilinear interpolation has an advantage over z interpola- 
tion in that it produces smooth fields. However, bilinear in- 
terpolation presents a number of problems: 

- The fields it produces will have a significantly different 
total mass than the fields produced by z interpolation. 
Our experience with monthly SMB fields over Green- 
land indicates that storms tend to leave large amounts 
of snow in a few localized areas. Bilinear interpolation 
tends to reduce the total amount of snowfall in these 
cases, causing potentially significant differences in the 
GCM model run. 

Variations in total mass caused by the choice of inter- 
polation procedure will not cause conservation prob- 
lems: we make this clear in Sect. 7. However, an in- 
terpolation scheme that produces a significantly differ- 
ent mass from the apparent “intent” of the GCM could 
cause biases. In particular, a scheme that makes snow- 
storms smaller than the GCM “intended” would pro- 
duce a negative bias in the equilibrium extent of the 
ice sheet. 

- It introduces significant numerical diffusion into the 
system. 

- The A ^ E transformation derived with bilinear inter- 
polation can introduce nonphysical artifacts (Sect. 12). 



Fig. 8. Setup for bilinear interpolation. The value in an iee grid eell 
(square) will be the sum of the values at the four nearest atmosphere 
grid eell eenters, weighted by A-longitude and A-latitude along the 
axes. 

- It is also not always clear how to extend bilinear in- 
terpolation to the case of non-rectangular atmosphere 
grids. This problem can also affect non-regular points 
in mostly regular grids (e.g., at the poles of a latitude- 
longitude grid). 

6.3 Elevation class interpolation 

One final choice for regridding is to define elevation classes 
as they were originally formulated (Leung and Ghan, 1998). 
In this approach, each atmosphere grid cell is grouped into 
sub-regions based on elevation bands. For example, an at- 
mosphere grid cell near the coast spanning elevations of 
0 m-600 m might be grouped into sub-regions of 0 m-200 m, 
200m-400m and 400 m-600 m. The grid E can then be 
thought of as having irregularly shaped grid cells formed by 
the intersection of atmosphere grid cell boundaries and ele- 
vation contours (Fig. 9). Area-weighted remapping is used to 
interpolate from the elevation grid to the atmosphere grid. 

Elevation class interpolation is equivalent to a specialized 
form of z interpolation, in which the I-D function fi{z) (for 
grid cell At) is interpolated as a discontinuous piecewise 
constant function (Fig. 10), rather than a piecewise poly- 
nomial (Fig. 4). Unless one expects significant discontinu- 
ities in topography or the elevation-flux relationship, such a 
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Fig. 9. Traditional elevation elass sehemes are equivalent to run- 
ning the iee surfaee model on an LO grid, where grid eell outlines 
are ereated by the interseetion of the atmosphere grid and elevation 
eontours. One sueh grid is shown in this figure. Note that grid eells 
extend only as far as the iee sheet; the grey line shows the Greenland 
eoast. 


zero-order interpolation scheme would be expected to give a 
less accurate version of than a higher order scheme such 
as z interpolation. Since we expect the elevation-flux rela- 
tionship to be smooth with elevation, we recommend the use 
of z interpolation over elevation class interpolation. 

6.4 Interpolating to LI grids 

We have outlined methods for the interpolation ^ G, as 
long as G uses an LO parameterization. The procedures need 
to be modifled for interpolation or ice grids using LI param- 
eterization. In a flnite element mesh (“grid”), fleld values are 
determined at mesh vertices and linearly interpolated within 
each triangular element (Zienkiewicz et ah, 2013). Any of 
the interpolation methods above may be used to determine 
the value of (x , y) at each vertex. Once vertex values have 
been interpolated on the vertices of a flnite element mesh, the 
value of f^(x, y) at all other points is fully determined. 


Surface Mass Balance 




• 


Atmospheric Mean 



^^0 1000 2000 3000 

Elevation (m) 


Fig. 10. Interpolated SMB funetion within one atmosphere grid eell 
when using elevation elasses. The resulting pieeewise eonstant in- 
terpolation is almost never preferable to the pieeewise linear inter- 
polation in Fig. 4. Traditional elevation elass sehemes are diseour- 
aged beeause they offer no benefit over first-order z interpolation. 


7 Transformation: elevation to atmosphere grid 

As shown in Fig. 2, a regrid step from the elevation to the 
atmosphere grid is required on every GCM time step. Once 
the modeler has chosen the transformation E G, we show 
here how to derive a transformation iox E ^ A that is con- 
sistent with it. 

To derive this transformation, consider one GCM time 
step, during which a flux fleld between the atmosphere 
and ice surface is computed on the elevation grid. That 
fleld will be regridded both to the interpolation and atmo- 
sphere grids. Conservation requires that the resulting fields 
are equivalent on A, i.e, 

E- ( 6 ) 

If the transformation E ^ G derived in Sect. 6 is repre- 
sented by the matrix M, then we can write f^= M/^ . Sub- 
stituting into Eq. (6), we get the requirement 

/A ^A (7) 

Equivalence on A can be tested by conservatively remap- 
ping both sides to A and then comparing (in this case, the 
left-hand side is already on A, so no remapping is required 
there). Denoting R as the matrix that conservatively regrids 
from G to A, we have the requirement 

/^=RM/^. (8) 

If we treat this as a definition for the transformation E 
A, we have derived a transformation for £ ^ A that is con- 
sistent with the transformation we already chose for ^ G. 

In other words (Fig. 11): we will regrid from £ to A 
by first regridding from to G (represented by the ma- 
trix M) and then use area-weighted remapping from G to 
A (represented by the matrix R). The GCM, which needs to 
use the transformation E ^ A, does not need to know any- 
thing about the ice grid used to derive that transformation. It 
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Fig. 11. The transformation from elevation grid to atmosphere grid is derived by first interpolating to the iee grid (represented by matrix 
M) and then using a eonservative area-weighted remapping step to the atmosphere grid (matrix R). The two steps may be eombined by 
eomputing the matrix produet RM. This eonstruetion ensures that a quantity eomputed on the elevation grid will have the same total mass 
when regridded to the iee grid or atmosphere grid. The two transformations from the elevation grid to the iee and atmosphere grids are said 
to be consistent. 


just needs to know the final matrix RM, which can be pre- 
computed via matrix multiplication. 

We have derived a transformation iox E ^ A that by def- 
inition is consistent with a previously chosen transforma- 
tion for ^ G. Our approach differs from previous efforts, 
which would start out with E ^ A and try to find a transfor- 
mation for ^ G that is consistent with it. 

7.1 A basis for the elevation grid 

Every vector space has basis functions, including the eleva- 
tion grid E. Our choice for the transformation E ^ A con- 
strains the basis functions we use on E. To find these basis 
functions, we expand on the central idea in the previous sec- 
tion, i.e., that we can determine properties of a field on the 
elevation grid by regridding to the ice grid. In this case, we 
define /^(x, y) to be equal to /^(x, y), where = M/^. 
In other words, we can evaluate /^(x, y) at a point by in- 
terpolating to G and then evaluating f^(x,y). This def- 
inition is consistent with our method for in the previous 
section. 

We now use this principle to obtain a formula for the basis 
functions of E. Substituting into Eq. (1), we get 


f^{x,y) = f^{x,y)=Mf^-g{x,y). (9) 

Rewriting in indicial notation, using the commutative 
property of multiplication, and swapping index letters, we 
get 

f{x,y) = f!^Mjigj{x,y). (10) 

We can use this, along with Eq. (1), to determine the basis 
functions for the elevation grid: 

6i{x,y) = Mjigj{x,y) (II) 

or, switching back to vector notation, 

e{x,y) = MJ g{x,y). (12) 


Equivalently, the basis function ei{x,y) is the function 
obtained if we set the /-^ = 1 and all other components of 
= 0, regrid to the interpolation grid, and then examine 
the resulting function f^(x, y). 

In general, these basis functions will be not orthogonal. 
Their exact form depends on choices made in choosing M, 
i.e., vertical and horizontal interpolation choices, as well as 
grid geometry issues. We have plotted some example basis 
functions in Fig. 12. 

Note that if elevation class interpolation is used (Sect. 6.3) 
and G is the exchange grid, then our scheme reduces to a 
traditional elevation class scheme, and basis functions will 
represent constant- value sub-grid tiles, which are orthogonal. 

7.2 Forward transformations: summary 

We have now defined three of the five regridding transforma- 
tions required by Fig. 2: E ^ G, G ^ A and E ^ A. We 
call these the forward transformations because they are lin- 
ear and can be represented by matrices. See Fig. 13 for a 
diagram of how these transformations may be used to regrid 
fields. 

We defined these three transformations in a consistent 
manner. We began by allowing the user choice in construct- 
ing the transformation from the elevation to interpolation 
grids, represented by matrix M: this makes sense because 
users desire choice for E ^ G. We then noted that standard 
area-weighted remapping, represented by matrix R, may be 
used to regrid from the interpolation grid to the atmosphere 
grid. From these two choices, transitivity requires that the 
transformation from the elevation to atmosphere grid is rep- 
resented by the matrix RM. Conservation is maintained be- 
cause E ^ G and E ^ A are consistent with each other, 
producing and f^ that are equivalent on every grid cell 
in A. 
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Elevation Basis Functions 


1150m 



Fig. 12. Unitless basis functions for the elevation grid E, constructed using 20 elevation points and z interpolation (the exchange grid was 
used as the interpolation grid). The grey box represents one atmosphere grid cell on the west coast of Greenland, with the coastline shown as 
black lines. White lines in the atmosphere grid cell represent elevation contours corresponding to each elevation point. The basis functions 
corresponding to elevation points at 950 m, 1 150 m and 1350 m are shown. Note that basis functions overlap and are not orthogonal. Because 
of the z interpolation, each basis function has maximum value at its corresponding elevation, but it has a non-zero support up to one elevation 
point away. 


Our choices for these transformations lead directly to a 
well-defined set of basis functions for the grid E, examples 
of which we plotted. Fields on this grid will be represented 
in terms of these basis functions - for example, surface- 
atmosphere fluxes and ice surface model state. 

8 Reversing the transformations 

We have defined three transformations so far, but Fig. 2 indi- 
cates that five are necessary for full functioning of the cou- 
pled system. We still need to derive appropriate procedures 
for the “reverse transformations” E and G ^ E, indi- 
cated in Fig. 13 as dotted lines. Because of the differences in 
dimensionality between the three grids (Sect. 2.3), M and 
RM are not invertible in a simple linear algebraic sense: 
A^Eis underdetermined and G ^ E is overdetermined. 
We must still And ways to compute these transformations that 
retain conservation over atmosphere grid cells. 

8.1 Transformation: atmosphere to elevation grid 

Suppose we have a flux field on the atmosphere grid A. 
We wish to regrid it to an “equivalent” flux field on the 
elevation grid E. In this case, / might represent precipitation 
or downwelling radiation from the atmosphere. 

The problem is underdetermined. Any solution for which 
fA j-E pg conservative. This is the same as requir- 
ing that 

RM/^ = /^. (13) 


Because of the many-to-one relationship between eleva- 
tion and atmosphere grid cells, our intuition tells us that 
may be constructed simply by repeating values of 
within each atmosphere grid cell. This is physically self- 
consistent. More precisely, if Ej e At, then we would like 
to set = /.^ . We define a simple linear transformation A 
that does this: /^ = A/^. 

Surprisingly, this definition for A ^ E is only conserva- 
tive in some cases. If RM computes each /-^ as a weighted 
sum only of values where Ej g Ai , then we say that RM 
is a local transformation - it uses only data from “within” 
an atmosphere grid cell to compute a value on that grid cell. 
In that case, it is easy to prove conservation by showing that 
RMA (see Eq. 13). The weights involved in RM 
do not matter, as long as they sum to 1 for each component 
of f^. 

Therefore, if RM is local, we can simply use A for our 
transformation A ^ E. Even better, it is easy to show that 
there will be no numerical dispersion in the round-trip trans- 
formation E ^ A ^ E because ARM is the identity matrix. 
The lack of dispersion on this round-trip is important because 
it is computed every atmosphere time step. 

8.2 Nonlocal RM 

If RM is not local, then in general A and will not be 
equivalent on A or even the entire ice sheet. We are forced to 
trade off between the most physically self-consistent value 
for = vs. something that conserves. We can use A, 
along with quadratic optimization, to guide us to such a com- 
promise. 
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We seek a vector such that 
RM/^ = /^. (14) 

We also wish that vector to be as close as possible to our 
intuition above. That is, we wish to minimize the quantity 

ll/^-A/^||2, (15) 

where our L 2 norm is weighted by the weight of each basis 
function. That is, we use a weight vector w where 

Wi = j ei{x,y)dA. (16) 

00 

This is a sparse quadratic optimization problem with 
equality constraints. A number of numerical packages can 
solve it; we used GALAHAD (Gould et al., 2003). In our 
tests, solution typically requires a fraction of a second on a 
single core. This is so fast that we have found no need to 
consider suboptimal solutions to this problem. 

This procedure introduces some numerical dispersion into 
the fully coupled system from nonlocal regridding opera- 
tions; by numerical dispersion, we mean movement of mass 
between adjacent atmosphere grid cells, thereby violating 
our desired property of conservation with each atmosphere 
grid cell. Considering the round-trip transformation E 
A ^ it is clear that this transformation is not local, both 
because RM is not local, and because the quadratic program 
set up for A ^ will be nonlocal. 

8.3 Practical issues for nonlocal RM 

GCMs are well equipped to deal with sub-grid tiles, each one 
occupying a fraction of the overall grid cell. The GCM will 
typically implement E ^ Ahy computing a weighted sum 
over sub-grid tiles in each grid cell, with weights based on 
each tile’s fractional area. This capability is used to imple- 
ment traditional elevation class schemes. 

If the RM matrix is local, then by definition the value 
on each atmosphere grid cell is a weighted sum of the ele- 
vation points in that grid cell. This is compatible with ex- 
isting GCM practice that assumes sub-grid tiles. To use the 
methods in this paper, the computation of £ ^ A inside the 
GCM does not need to be replaced, it only needs to be fed 
a new set of weights. Even though elevation points do not 
have well-defined areas, the weight for elevation point j can 
be thought of as the “fractional area” that an elevation point 
contributes to its containing atmosphere grid cell. The GCM 
can be coded as if traditional elevation classes were being 
used, even if the user has chosen a more numerically accu- 
rate form of vertical interpolation in the construction of the 
M matrix. 

Things get more complicated for the GCM if the RM ma- 
trix is not local. Instead of using a simple set of weights, 
the GCM will have to be multiplied by a general sparse ma- 
trix RM when computing E ^ A. An MPI (Message Pass- 
ing Interface) gather will be required to compute A ^ E. 





Fig. 13. The five grids used in the eoupling problem, and transfor- 
mations between them. The (linear) interpolation step M, from E 
to G, is ehosen by the user. The (linear) transformation R, from G 
to A, is an area-weighted remapping step. These two fully eonstrain 
the transformation from F to A as RM. If I is LO-parameterized, 
then G is the exehange grid between A and /, and X and X' are area- 
weighted remapping transformations. If I is LI -parameterized, then 
G is equivalent to /, making X and X' the identity. P is a diagonal 
resealing operation between A' and A. Reverse transformations re- 
quired by the eoupled system are shown as dotted lines. 

Finally, some GCMs use implicit schemes for ice surface- 
atmosphere coupling (Best et al., 2004), requiring a matrix 
inversion on every time step. If RM is nonlocal, then a large, 
sparse, global matrix inversion will be required, rather than a 
number of small, local inversions. 

Because of the significant practical complications that 
arise with the use of a nonlocal RM matrix, most users find 
it simplest to use a local RM matrix if at all possible. 

8.4 When is RM not local? 

The use of a nonlocal RM transformation can be an added 
burden to the GCM developer. It is therefore important to 
delineate cases in which RM is not local. 

For LI ice grids, as used with a finite element ice flow 
model, RM will be nonlocal: ice elements that straddle two 
atmosphere grid cells will by necessity involve elevation 
point values from two atmosphere grid cells. This will make 
R slightly nonlocal. 

Even for LO grids, M could be nonlocal, depending on the 
choice of interpolation schemes for E ^ G. Bilinear inter- 
polation is inherently nonlocal, whereas the other interpola- 
tion strategies mentioned above are local. 

For LO ice grids with a local interpolation scheme, RM can 
be made local, as long as G was chosen to be the exchange 
grid (see Sect. 4.4). In this case, M will be local - because 
each exchange grid cell attains its value from elevation points 
in just one atmosphere grid cell. R will be local as well, for 
the same reason. 
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8.5 Transformation: ice/interpolation to elevation grid 

Suppose we have a flux fleld on the interpolation grid 
G. We wish to regrid it to an “equivalent” flux fleld on 
the elevation grid E. As before, we will set up a quadratic 
optimization problem. 

It would be nice if we could And such that M/^ = 
/^. This will not usually be possible, since G has far more 
degrees of freedom than E. Instead, we will minimize the 
quantity 

l|M/^-/«||2 (17) 

while we maintain conservation on each atmosphere grid cell 
(where R deflnes area-weighted remapping from G to A; see 
Fig. 13): 

RM/^ = R/‘^. (18) 

This quadratic optimization problem may be solved with 
the same methods as in Sect. 8.2. We weight components of 
our L 2 norm by the integral of each basis function in G. 

9 Exchange grid or ice grid? 

So far, we have deflned our transformations in terms of the 
interpolation grid G - which could be either the ice grid or 
the exchange grid. But real ice flow models operate on the ice 
grid, and transformations must ultimately transform to/from 
that grid. If we have chosen G as the exchange grid, we need 
to extend our transformations above to regrid to/from the ice 
grid. We do this by using area-weighted remapping as neces- 
sary between G and I. 

More formally, we construct transformations for X : G ^ 
I and X' : I ^ G using area-weighted remapping. We can 
then represent E ^ I as the matrix product XM. Similarly, 
we can construct I ^ E hy first computing f^=X!f^ and 
then using the reverse transformation G ^ E to obtain 
(see Fig. 13). 

Note that the transformation X : G ^ / is not conservative 
on A. For this reason, the transformation E ^ 1 represented 
by XM is not conservative on A either, although it is con- 
servative overall: “mass” lost from one atmosphere grid cell 
will be gained by neighboring cells. This will cause numer- 
ical dispersion when these transformations are used in the 
fully coupled system (Fig. 2). 

Assuming the ice grid is LO parameterized, we will now 
address the practical issue faced by the user, namely whether 
to choose the exchange grid or ice grid as G. If one uses the 
exchange grid for G in conjunction with a local vertical inter- 
polation scheme, then the RM matrix will be local. This has 
a number of advantages. It eases implementation (Sect. 8.3). 
There will be no numerical dispersion for the transformation 
E ^ A, used every atmosphere time step. The reverse trans- 
formation A ^ E, also used every atmosphere time step. 


will be trivial. However, there will be numerical dispersion 
for the transformation E ^ I, which is used once every ice 
coupling time step. 

If one directly uses the ice grid for G, then the RM ma- 
trix will not be local, producing numerical dispersion for the 
transformation E ^ A and complicating the reverse trans- 
formation A ^ E. However, there will be no numerical dis- 
persion for the transformation E ^ I. 

On balance, we recommend keeping the RM matrix local 
if possible. Not only does this simplify implementation, it 
moves unavoidable numerical dispersion away from the at- 
mosphere time step to the less frequent ice coupling time 
step. 

9.1 Exchange grids for finite elements 

The discussion in Sect. 9 assumes an LO ice grid. The results 
can be extended to LI or higher order parameterizations: fi- 
nite element ice models, for example. In this case, the con- 
cept of exchange grid is not meaningful, since exchange grids 
are by definition LO. 

Instead of the exchange grid, one can choose G to be just 
about any LO grid, of resolution at least as high as the ice 
mesh, where grid cells do not cross atmosphere cell bound- 
aries: we will call this a generalized exchange grid. The 
user would then need to develop an appropriate conserva- 
tive transformation for G ^ /, whereas the transformation 
for / ^ G would remain an area-weighted remapping (see 
Appendix C). Details of this scheme for a particular ice mesh 
parameterization are left to the reader. 

As in Sect. 9, the transformation G ^ / would cause some 
numerical dispersion, due to the geometry of the ice mesh. 
Additional dispersion would be introduced because of mis- 
match between the LO grid G and the higher order mesh I 
- although this could be reduced by making G finer. The 
trade-offs of choosing G to be a generalized exchange grid 
vs. the ice grid are the same as in Sect. 9: numerical disper- 
sion is introduced in the ice coupling time step, in exchange 
for simplifying implementation and eliminating dispersion in 
the atmosphere time step. 


10 Regridding examples: local RM 

Having shown how to compute all five required regridding 
transformations, we now demonstrate them working within 
a realistic GCM context. We set up a test using the GISS 
2°x2i° atmosphere grid (Schmidt et al., 2006), overlap- 
ping with the SeaRISE 5 km LO grid (Bindschadler et al., 
2013) and 40 elevation points, spaced every 100 m from Om 
to 4000 m. We chose the exchange grid as G (see Appendix A 
and Fig. 3 for notation conventions). We ran GISS ModelE in 
this configuration, using fixed sea surface temperatures cor- 
responding to the years 1996-2005. The ice surface model 


Geosci. Model Dev., 7, 883-907, 2014 


www.geosci-model-dev.net/7/883/2014/ 


R. Fischer et al.: Conservative regridding for ice-atmosphere coupling 


897 



Fig. 14. July SMB computed by ModelE on the elevation grid E (a), and regridded to the atmosphere grid A via the exchange grid (b) or to 
the ice grid I using z interpolation (c). For plotting purposes, each elevation point in (a) is assigned to a nearby region of similar elevation; 
the value of the elevation point is then plotted in its corresponding region. Bold numbers on the color scale indicate the extreme values of the 
plot. Elevation contours are plotted. 


was run on the elevation grid, producing monthly averages 
of SMB over Greenland, which we label (Fig. 14a). 

As expected, the broad pattern shows strong melting in 
narrow bands at low elevations, along with weak accumu- 
lation at high elevations. Atmosphere grid cell boundaries 
are prominent because precipitation - the primary source 
of positive SMB - is not downscaled to sub-grid resolu- 
tion (Leung and Ghan, 1998). The small oscillations at high 
elevation are artifacts introduced by ModelE ’s ice surface 
model. 

10.1 Example: elevation to atmosphere grid 

Figure 14b shows the original field computed on the 
elevation grid and regridded to the atmosphere grid A, via 
the transformation =RM/^. Note how low-elevation 
ablation regions become broader and weaker (compared to 
Fig. 14a), whereas the interior of the ice sheet remains about 
the same. This demonstrates the benefits obtained through 
the use of elevation points, as compared with running the ice 
surface model on the atmosphere grid (van den Broeke et ah, 
2008). 

By definition, the transformation RM is conservative on 
A: we can evaluate conservation properties of other transfor- 
mations by comparing to Fig. 14b. 

10.2 Example: elevation to ice grid 

Figure 14c shows the original field regridded to /, via 
f^= XM using z interpolation for M. This plot looks 


like a smoothed version of Fig. 14a. Atmosphere grid cell 
boundaries are still visible because z interpolation is local. 

The transformation X : G ^ / (Fig. 13) introduces a small 
amount of nonlocality. And although it is conservative over 
the ice sheet in general, it is not conservative over A. This 
can be seen (Fig. 15) by regridding to A and comparing 
with the results of Fig. 14b. This plot quantifies the amount 
of numerical dispersion the simulation will encounter every 
month when preparing SMB input for the ice model. This 
dispersion is caused by ice grid cells that overlap more than 
one atmosphere grid cell. 

In most areas, numerical dispersion is low, less than 1 %. 
That is because most atmosphere grid cells are overlapped by 
many ice grid cells, with only a few lying on an atmosphere 
grid cell boundary. However, numerical dispersion can be 
significant for atmosphere grid cells that just nick the edge 
of the ice sheet, where a high proportion of their overlapping 
ice grid cells also overlap with a neighboring atmosphere grid 
cell. 

The numerical dispersion encountered here is inherent in 
the regridding problem itself, rather than our approach to that 
problem. Short of using an ice grid whose grid cells do not 
overlap atmosphere grid cells, we see no obvious way to 
eliminate this issue. However, we do not believe it to be a 
serious problem: the total area of ice sheet affected by it will 
be small. 
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Fig. 15. Difference in July SMB obtained with E ^ A, using two 
choices for the interpolation grid G: either the exchange grid or the 
ice grid. Differences (exchange minus ice) are due to ice grid cells 
that overlap more than one atmosphere grid cell. 


It is important to point out that although X' is not (quite) 
conservative over A, it is still conservative over the ice sheet 
in general. We have verified this numerically in our exam- 
ples. 

10.3 Example: ice to elevation grid 

Since we have not yet developed the surface boundary con- 
dition described in Sect. 2.1, we do not have realistic fields 
on I to try regridding to E. We will therefore test I ^ E at 
this point using f^= XM as input. 

We test the transformation in two steps. First, we test 
E using = M/^ as input. Then, we test 1 ^ E 
using f^= XM as input. This allows us to evaluate the 
G ^ E transformation separately from confounding disper- 
sion factors involved in G ^ I . 

We computed = [G ^ £](/^) where = M/^. 
Theoretically, f'^ and should be equal. The conservation 
property imposed by the quadratic program held: we found 
that f'^ to machine precision. However, we 

found measurable differences between f'^ and (Fig. 16). 
Although differences are small in most cases, one area has a 
difference of more than 8 out of 3 7 mm day resulting in 
no more than one digit of precision. Differences tend to be 



■ I I 
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Fig. 16. Example of numerical dispersion caused by an ill-posed 
quadratic optimization problem. July SMB on elevation grid E was 
regridded to the exchange grid. The reverse transformation was then 
used to recover the original SMB on E. Plot shows the difference 
between the result and the original. In theory, the two should be 
exactly the same. Although mass is conserved, differences appear 
because the associated quadratic optimization problem is underde- 
termined in some areas. 


greatest in sparsely populated atmosphere grid cells on the 
edge of the ice sheet. We conclude that the answer to the 
quadratic program posed in Sect. 8.5 is poorly constrained. 
A different numerical solver than the one we used might in 
theory result in a better match between f'^ and /^. 

But this is not a problem in practice: the goal of G ^ 
is to obtain a physically plausible field on E with the correct 
conservation properties. The preliminary tests presented in 
this section are consistent with that goal. 

The transformation / ^ £ is constructed by transforming 
I ^ G ^ E (Sect. 9). We used this to compute = [/ ^ 
where = XM /^. The conservation property im- 
posed by the quadratic program held: we found that f'^ was 
equivalent to on A to machine precision. 

However, differences between f'^ and are even 
greater in this case: compare Fig. 17 to Fig. 16. This is be- 
cause the transformation X : G ^ / is not conservative on 
A. The end result of transforming E^G^l^G^E 
will produce an that is significantly different from the 
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Fig. 17. Example of numerical dispersion in the transformation 

I ^ E, caused by ice grid cells that overlap more than one atmo- 
sphere grid cell. July SMB on elevation grid E was regridded to the 
ice grid /. The reverse transformation was then used to recover the 
original SMB on E. Plot shows the difference between the result 
and the original. Although mass is conserved, differences are due 
mainly to ice grid cells that overlap more than one atmosphere grid 
cell, making it impossible to recover the exact original SMB. The 
same color scale is used as in Fig. 16. 

original . Differences are greatest for elevation points near 
an atmosphere grid cell boundary. 

II Regridding examples: almost local RM matrix 

In Sects. 8.3 and 9, we discussed the trade-offs between using 
the ice vs. exchange grid for the interpolation grid G. Having 
shown in Sect. 10 an example of our transformations using 
the exchange grid as G, we now show how things change if 
the ice grid is chosen for G by re-running the above exam- 
ples. In this case, the RM matrix will be mostly local, except 
for nonlocality caused by ice grid cells that intersect more 
than one atmosphere grid cell. We say that RM is almost lo- 
cal. 

In general, differences in results here vs. those in Sect. 10 
are insignificant. However, the choice of I as interpolation 
grid does change the basis functions used for E, yielding a 
system in which the transformation £ ^ / is now conserva- 
tive over A. 


With a change in £ ^ the transformation E ^ A is 
also changed. We calculated the difference between com- 
puted using this method vs. computed in Fig. 14b: this 
difference is exactly the same as Fig. 15. In practice, this dif- 
ference does not present a real problem: it is impossible to 
say which version of is more “correct.” Both maintain 
conservation over the ice sheet. As predicted in Sect. 9, grid- 
related dispersion in I ^ E is eliminated. Thus, the grid- 
related “errors” demonstrated in Fig. 17 are eliminated. 

Finally, the use of a nonlocal RM matrix requires use of 
the algorithm described in Sect. 8.1 for A^E. Fig. 18a 
shows a July precipitation field produced by ModelE, for 
the same month as the SMB fields above. ModelE currently 
does not downscale precipitation to sub-grid resolution (Le- 
ung and Ghan, 1998): it therefore yields a precipitation field 
on A. 

We then computed = [E ^ A](p^), which could 
serve as input to the land surface model. Recall that Ap^ 
is the most physically self-consistent value for p^ , but that 
we choose a different value for p^ in order to maintain con- 
servation. 

Figure. 18b shows the difference between the two, p^ — 
Ap^. Differences are typically in the range [—0.1,0.!] 
(about 2 %). They tend to be relatively constant within each 
atmosphere grid cell and show no discernible pattern be- 
tween grid cells. In particular, differences are not related 
to elevation. Any errors introduced by this scheme will be 
dwarfed by other precipitation errors in the model. 

In this example, we used an LO ice grid to generate a pro- 
totypical almost-local RM matrix. We expect similar results 
when using an El ice grid because the nonlocality in the LO 
case was caused by a small number of ice grid cells that over- 
lap more than one atmosphere grid cell. A similar situation 
exists with any type of higher order mesh. 

We conclude by considering, in the case of an LO ice grid, 
whether / or G is a better choice for the interpolation grid. In 
many cases, the use of a nonlocal RM matrix requires signif- 
icantly more effort in GCM model development. However, 
our experience shows fewer “surprises” in the transforma- 
tions when using I as the interpolation grid. In the end, we 
expect either choice to yield serviceable results that conserve 
mass and energy in the regridding. 

12 Regridding examples: nonlocal RM matrix 

In the above sections, we considered how the transformations 
in this paper work when the RM matrix is local or almost lo- 
cal. Here, we consider the case in which RM is significantly 
nonlocal - for example, if bilinear interpolation were used 
for the transformation E ^ 1. We would expect the “correc- 
tions” and dispersive properties that appear in the case of an 
almost-local RM matrix to be significantly larger. 
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Fig. 18. Nonphysical changes in precipitation field introduced by A ^ F in order to ensure conservation. Panel (a) shows a July precipitation 
field computed by ModelE on the atmosphere grid. This field was regridded to the elevation grid E using an almost-local RM matrix for 
E ^ A. It produced a slightly different result than the physically intuitive procedure of “repeating” precipitation values for elevation points 
in each atmosphere grid cell — which is not conservative. Those differences are plotted in panel (b). 


We encountered numerous problems when we attempted 
the use of a significantly nonlocal RM matrix, constructed 
using bilinear interpolation for E ^ I. Most significantly, 
the A ^ F transformation produced elevation points of neg- 
ative precipitation when regridding a typical precipitation 
field over Greenland (Fig. 19) - an artifact that we believe 
would be unacceptable to the majority of modelers. It might 
be possible to find a solution to this problems. In the mean- 
time, it is simpler to avoid nonlocal interpolation schemes 
such as bilinear interpolation. 

13 Independent elevation grid 

The elevation grid E is constructed by adding elevation 
points to each grid cell of an underlying horizontal grid, 
which we will denote F^. So far we have assumed that 
F^ = A, meaning the elevation grid E is derived from the 
atmosphere grid A. Some GCMs use a horizontal layout for 
E that is not related to A. In this section, we extend our meth- 
ods to address that case. 

The first problem is a choice of the interpolation grid G. 
The idea of locality in RM no longer makes sense when 
E^ ^ A. For that reason, G = / is the right choice, X and 
X' will both become identity transformations. The transfor- 
mation E ^ I does not need to change, nor does / ^ A: 


none of these rely on any specific relationship between F^ 
and A. Similarly, F ^ A = RM can be computed as before. 

The only other thing that must change is the construction 
of A ^ F: the intuitive definition for A (Sect. 8.1) no longer 
makes sense. Instead, we construct an intuitive regridding 
operation F as follows. Given /^, first regrid to using 
area- weighted remapping. Then convert to using a 

“repeat” operator akin to A above. We can now apply the 
quadratic programming-based regrid operator developed in 
Sect. 8.2, using F instead of A in Eq. (15). 

The use of F^ / A introduces numerical dispersion into 
the system. This can be seen by evaluating the locality of the 
round-trip transformation E ^ A ^ E. This numerical dis- 
persion is a fundamental consequence of the use of a different 
underlying grid for the ice surface and atmosphere models. 

14 Regridding in elevation space 

We have tacitly assumed so far that there is one single fixed 
elevation grid E with one fixed set of basis functions. This is 
not the case in a changing climate, since the basis functions 
used for E depend on ice elevation and extent (Sect. 7.1). 
The user might wish to explicitly move elevation points as 
well - for example, to track a mountain glacier as it moves up 
or down. When the vector space E changes, the ice surface 
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Fig. 19. Implausible precipitation field produced hy E when 
bilinear interpolation is used. July precipitation (Fig. 18a) was re- 
gridded to the elevation grid E using the A ^ E reverse transfor- 
mation, where the matrix for F ^ A is constructed using (nonlo- 
cal) bilinear interpolation. Note the unphysical artifacts (negative 
precipitation) for some elevation points. 


model state must be regridded from the old elevation grid to 
the new elevation grid. We address that issue in this section. 

Assume two elevation grids, an “old” grid E and a “new” 
grid F. We wish to conservatively regrid a field to /^. In 
this case, / will not be a fiux variable, but rather a conserved 
state variable of the ice surface model: snow depth, water 
fraction, etc. 

This problem can be approached by examining the system 
of grids and transformations available in Fig. 13 when one 
has multiple elevation grids (see Fig. 20). By “following the 
arrows,” the most direct way to transform to is to first 
regrid to the interpolation grid, computing = 

Then use the procedure in Sect. 8.5 to compute /^. Because 
all these transformations are conservative on the atmosphere 
grid A, the end result will be equivalent to on A. 

14.1 Conserved model state 

In previous sections, / represented fiuxes between models, 
which are conserved. In this section, / is a model state vari- 
able of the ice surface model. In order for this regridding 



Fig. 20. When ice extent, ice topography or elevation points 
change, the basis functions for the elevation grid E change along 
with it. Ice surface model state, which exists on E, must be regrid- 
ded to the new set of basis functions. Shown is a grid system that 
can serve as a map for this regridding: E is the old elevation grid, 
E is the new one and G is the interpolation grid (same for old and 
new). Ice surface model state may be regridded from F to F by first 
regridding E ^ G, then G ^ E . Note that this diagram is a simpli- 
fied version of Fig. 13 in which two different elevation grids have 
been accounted for. 

procedure to be physically meaningful, the model state must 
be expressed in terms of conserved quantities. 

Not all ice surface models are formulated in terms of con- 
served quantities. In this case, the regridding procedures may 
still be used, as long as the ice fiow model can be converted 
to/from a form that is expressed in conserved quantities. For 
example, an ice surface model might track the temperature, 
mass and water fraction of the top layer of ice. Temperature 
is not conserved, so this regridding procedure cannot be used 
directly. However, model state can be converted to enthalpy 
and mass alone (Aschwanden et al., 2012). These quantities 
are conserved and can be correctly regridded with conserva- 
tive transformations. After regridding to the new elevation 
grid, model state can then be converted back to the original 
nonconserved parameterization. 

14.2 When to regrid 

Regridding in elevation space is not just required if one 
changes elevation points, but also any time that the transfor- 
mation M for F ^ / changes. Any time that happens, the 
regridding described in Sect. 14 should be used. Relevant 
events that change M include 

1. changes in ice topography: for example, as an ice sheet 
infiates or defiates due to changing climate; when an 
ice grid cell changes elevation, the weights in M used 
to compute it will change; this change will either be 
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Fig. 21. The Glint2 workflow used to eompute the regridding opera- 
tions required by the fully eoupled system (Fig. 2). Glint2 produees 
regridding operations based on a variety of faetors: atmosphere and 
iee grid geometry, iee topography and extent, and levels ehosen for 
the elevation grid. Glint2 must reeompute the operations when any 
of these faetors ehanges. 

continuous (as in z interpolation) or will jump at cer- 
tain thresholds (as in elevation class interpolation). 

2. changes in ice sheet extent, as an ice sheet grows or 
melts; for example, if an ice sheet shrinks, then some 
ice grid cells will no longer participate in the regrid- 
ding process, and their associated coefficients in M 
will turn to zero. 

3. changes in the elevation or number of elevation points. 

15 Glint2 coupling library 

Here we describe Glint2, an open-source implementation 
of the transformations developed in this paper. Glint2 is a 
GCM-ice flow model coupling library whose core ffinction 
is to compute the flve transformations required for a coupled 
GCM-ice flow model system (Fig. 2). These transformations 
are computed based on a variety of factors: atmosphere and 
ice grid geometry, ice topography and extent, and elevation 
levels chosen for the elevation grid (Fig. 21). 

Glint2 does not just compute these transformations and 
provide them as a library, it also provides an application pro- 
gramming interface (API) for GCMs to use to couple with 
ice flow models. As a realization of the mediator design pat- 
tern (Gamma et ah, 1994), Glint2 is able to shield the GCM 
from a number of details of the coupling. Although it uses a 
different codebase, Glint2 has the same purpose as GLINT 
(Glimmer Interface; Rutt et al., 2009), namely to build a 
coupling library between GCM and ice models. The GCM 
programmer who wishes to couple their GCM with an ice 
model need only do the following: 

1 . implement an elevation points scheme, and move the 
ice surface model to the elevation grid. 

2. request the RM matrix from Glint2, and then apply 
it as needed during normal model run; this might be 
more complex than it sounds, depending on the GCM’s 


system of domain decomposition; storage and applica- 
tion of the RM matrix is simplifled if the GCM author 
knows in advance that it is local. 

3. ask Glint2 to regrid ice surface model inputs from 
A ^ as needed; this step is only required if RM 
is nonlocal; otherwise, the transformation A is simple 
enough for the GCM to do on its own. 

4. accumulate ice surface fluxes on the elevation grid, and 
pass them to Glint2 every coupling time step. 

5. apply flelds returned from the ice flow model via 
Glint2; these flelds are returned on the elevation and 
atmosphere grids as appropriate, eliminating the need 
for the GCM to regrid them. 

Note that the GCM does not need to know anything about 
the ice grid or ice flow model. Interface code is added to 
Glint2, not the GCM, for every ice flow model one wishes 
to support. Over time, we expect the number of supported 
GCMs and ice flow models to increase, according to re- 
searcher demand. This structure is useful because it will give 
practitioners a way to try different ice models with a GCM 
fairly easily. 

15.1 Adoption issues 

We expect that Glint2 could be useful to anyone with a 
GCM who wishes to couple it with an ice flow model. How- 
ever, many GCMs have centralized regridding strategies that 
Glint2 does not really fit into. We do not believe this should 
be a significant barrier to adoption for two reasons: 

- In most cases, the elevation grid will have unusual 
“customized” basis functions (Fig. 12). Centralized 
GCM regridding schemes are not typically equipped 
to regrid to/from the elevation grid. Nor are they 
equipped with the algorithms required for the reverse 
transformations. 

- Glint2 hides all details of the ice grid from the GCM, 
communicating with the GCM with flelds on the atmo- 
sphere and elevation grids. Because all ice grid-related 
issues are encapsulated in Glint2, it should not matter 
to the GCM how or even whether regridding to/from 
the ice grid is accomplished. As far as the GCM is con- 
cerned, the ice flow model might as well be running on 
the elevation grid. 

Another barrier to adoption is the fact that Glint2 is pack- 
aged as a library. Many GCM projects are reluctant to add ad- 
ditional library dependencies, due to the complications such 
dependencies introduce in the build process. We do not be- 
lieve this should be a serious problem because coupling with 
an ice model already involves the use of outside libraries. 
With or without Glint2, the GCM must still manage addi- 
tional external dependencies. 
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15.2 Availability 

Glint2 is written C++ and designed to couple with GCMs and 
ice flow models written in Fortran 90/2003, C or C++. It also 
comes with a Python interface, making it easy to test and plot 
sample regridding problems before incorporating a coupling 
strategy into a GCM. Glint2 source and documentation are 
available for download (Fischer, 2013). 

16 Discussion and conclusions 

This paper focuses on a system of conservative regridding 
strategies needed to support tight two-way coupling between 
a GCM and an ice flow model in the context of elevation- 
based downscaling in the GCM. Past efforts at one-way cou- 
pling have produced downscaling methods that provide SMB 
flelds that match remarkably well with observations and with 
regional models. However, these efforts used an inconsistent 
set of transformations, which would result in nonconserva- 
tion of mass and energy in a two-way coupled system. 

In order to achieve consistency, we began by recognizing 
the elevation grid as an integral part of the coupling problem, 
along with the atmosphere and ice grids that had previously 
been considered. We observed that the ice surface model runs 
on the elevation grid and that the downscaling step (from el- 
evation to ice grid) is actually another form of regridding. 
We have therefore transformed a coupling problem involv- 
ing two models and two grids into one with three models and 
three grids. 

We analyzed the regridding transformations needed in 
a typical two-way coupling and determined that live such 
transformations are required: three “forward” transforma- 
tions and two “reverse” transformations. We then set out to 
develop a consistent implementation for these flve transfor- 
mations, starting with the forward transformations. 

We observe that conservation of mass and energy requires 
consistency between the transformations from the elevation 
to ice grid {E ^ I) and elevation to atmosphere grid {E ^ 
A). We achieve consistency by allowing the user to choose 
E ^ I and then constructing a transformation for ^ A 
that is consistent with the user’s choice. This is a good ap- 
proach because it allows the user freedom in choosing a 
downscaling transformation for ^ /. Our transformations 
imply a set of basis functions for the elevation grid, which 
we demonstrated in plots. 


We then addressed the reverse transformations, using the 
notions of consistency developed in the previous sections. 
Problems of underdeterminism and overdeterminism, caused 
by mismatches in dimensionality between grids, are ad- 
dressed by using quadratic optimization for these transfor- 
mations. 

The result is a set of flve conservative transformations 
needed to support a two-way coupling of GCMs and ice 
models. We deflned a property of our E ^ A transformation 
called locality. We showed a number of theoretical and prac- 
tical benefits if that transformation is local. We implemented 
all flve transformations and demonstrated that they work in 
practice on realistic input flelds. 

Although three of those flve transformations are well- 
known in the literature (Ramshaw, 1985; Lipscomb et al., 
2013), the reverse transformations are new. Most impor- 
tantly, we showed how to choose a set of grids, basis func- 
tions and regridding schemes so that all five regridding op- 
erations are conservative. Note that the elevation grid uses 
a “custom” set of basis functions. This is a significant step 
beyond past efforts that have focused on conservative regrid- 
ding to LO pre-chosen grids (Ramshaw, 1985). 

Our downscaling transformation from the elevation to ice 
grids is taken from previous one-way studies (Lipscomb 
et al., 2013). Our results “look” almost identical, except for 
subtle differences in the other related transformations needed 
to ensure conservation. We therefore provide practitioners 
with a way to bring already proven downscaling techniques 
into a conservative two-way coupled setting. 

We went on to package these transformations in Glint2, 
a coupling library. As a piece of software, Glint2 serves as 
a mediator (Gamma et al., 1994) between the GCM and ice 
flow model, insulating each from the specific details of the 
other. This will simplify the coupling of GCMs with multiple 
ice models, as needed to support future research. 
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Appendix A 

Mathematical conventions 

We use the following mathematical conventions in this paper: 

1 . A “grid” or “parameterization” is a set of basis func- 
tions that may be linearly combined to produce func- 
tions over a 2-D domain. Grids are denoted by non- 
bold capital letters: A for the atmosphere grid, I for 
the ice grid, E for the elevation grid and G for an inter- 
polation grid. Basis functions are denoted by non-bold 
lower case letters: ai (x, y) is a basis function for A. At 
times, the set of basis functions may be represented as 
a bold- face vector: a(x , y) = [a\(x , y) , . . . , a^ix , y)]. 

2. Subscripts are used to indicate elements of a vector. If 
/ is a vector, then fi is the value of the /th index of /. 

3. If A is an LO-parameterized grid and A/ one of the grid 
cells, then | A/ 1 is the area of that grid cell. 

4. Vector values are indicated in boldface, scalars in non- 
bold. 

5. Vectors are used to construct a 2-D function (“field”) 

within the context of a particular grid. For example, 
suppose that the atmosphere grid A uses the basis func- 
tions «(x, y) = [< 2 i(x, y), . . . , y)], and we have 

an ^-dimensional vector with components /.^. 
That vector represents the function: 

f^(x,y) = f^-a{x,y). (Al) 

6. Since this paper is about regridding, we need to talk 
about the “same” field in different grids. Superscripts 
are used for this: if the vector denotes the vector 
representing a field on the ice grid, then denotes a 
vector representing the “same” field on the atmosphere 
grid. 

7. Similarly, f^(x,y) denotes the function implied by 

the vector and the vector of basis functions «(x, y). 

8. We use arrows to talk about transformations between 
grids, which transform fields on one grid to fields on 
another. Since fields on a grid are represented by vec- 
tors, the transformations are functions from one vector 
space to another. For example, E ^ A can be read as 
“the transformation from the elevation grid to the at- 
mosphere grid.” Many but not all of these transforma- 
tions are linear and can be represented by matrices. 

9. If a transformation is nonlinear, we use functional no- 

tation with square brackets to denote an application of 
that transformation: = [A^ £](/^). 

10. We use indicial notation in one place. 

11. See Fig. 3 for definition of the symbols used through- 
out the text. 


Appendix B 
Integration of LO grid 

Suppose we have a region B and an LO-parameterized grid 
G with grid cells Gi , . . . , G„. The basis function gi (x,y) is 
equal to 1 inside of grid cell Gi and zero elsewhere. We wish 
to compute 

I gi(x,y)dA = \BDGi\. (Bl) 

B 

This is simply equal to the area of intersection of B and 

G,. 

If we wish to regrid to an LO grid A, we must do the above 
computation multiple times, setting B to every grid cell Ay g 
A. This creates an overlap matrix L where L/y is equal to the 
overlap between source grid cell Gi and destination grid cell 

The overlap matrix is directly related to the exchange grid 
between G and A: every non-zero element of L/y is equal to 
the area of one grid cell in the exchange grid. The exchange 
grid may be computed as described in Appendix D. The Sur- 
veyor’s Formula (Braden, 1986), a special case of Green’s 
Theorem, can then be used to compute the overlap matrix 
from the exchange grid. 

Appendix C 
Integration of LI grid 

Appendix B shows how to integrate over areas on an LO ice 
grid. But if an LI finite element ice flow model is used, one 
will need to do this over an LI ice grid. Here, we show how to 
integrate LI finite element basis functions over an arbitrary 
area B. 

An LI finite element mesh is made up of triangular ele- 
ments, where the value of a function is defined at the vertices 
of the triangles. Values inside each element are interpolated 
based on the values at the vertices. 

Each basis function Ni (x , y ) in a finite element mesh cor- 
responds to a vertex (Zienkiewicz et al., 2013). This function 
is non-zero only in the triangular elements in which the ver- 
tex participates. We define a sub-basis function Nij(x, y) to 
be equal to Ni(x,y) within triangular element j, and zero 
elsewhere. Thus we define 

Ni{x,y) = J2Nij(x,y). (Cl) 

j 

Within its element, a sub-basis function takes the shape of 
a plane, having the functional form 

Nij (x,y) = a -\-bx cy. (C2) 

The coeflicients a, b and c are functions of the locations 
of the vertices, not of the values assigned to those vertices. 


Geosci. Model Dev., 7, 883-907, 2014 


www.geosci-model-dev.net/7/883/2014/ 


R. Fischer et al.: Conservative regridding for ice-atmosphere coupling 


905 


Using the techniques from Appendix D, one can compute the 
polygon Pj as the intersection between element j and region 
B. Green’s Theorem may then be used to compute 

j Nij{x,y)dA. (C3) 

Pi 

We can then apply our definition of sub-basis functions to 
achieve our goal: 

j Ni{x^y)dA = ?/ Nij{x,y)AA. (C4) 

B J ftj 

This section has provided just an outline of the process. 
The algebra can become complex at times, and a symbolic 
computation system such as Maxima can be useful. But in the 
end, integration of a vector over an area B is computed as 
a linear combination of the elements of /^, just as with LO 
grids. 

Appendix D 

Computing the exchange grid 

Appendices B and C can be applied repeatedly to compute 
conservative regridding matrices from an LO or LI grid to an 
LO destination grid. Both procedures assume a way to com- 
pute polygonal intersections between two sets of polygons - 
also known as an exchange grid (Balaji et al., 2006). 

This task is simple in principle, using modern computa- 
tional geometry packages such as CGAL (2013). If A and G 
are two sets of polygons, we explicitly construct each poly- 
gon in A and G and then compute pairwise intersections be- 
tween the two sets (Chin and Wang, 1983). This algorithm 
provides not just integration formulas, but also the actual 
polygonal outlines of the exchange grid. 


If an appropriately robust polygon intersection algorithm 
is used, our procedure can deal with nonconvex polygons and 
other possible irregularities. This is not just a theoretical is- 
sue: latitude-longitude grid cells commonly used in GCMs 
are not convex in spherical or planar geometry. In other cases, 
practitioners might wish to use grid cells consisting of mul- 
tiple disjoint polygons. 

With issues of polygonal intersection taken care of by pre- 
packaged algorithms, the main challenge here is to find those 
intersections in a scalable manner. The naive algorithm is to 
write a nested loop, requiring \A \ x |G| iterations. This algo- 
rithm, with 0(n^) complexity, takes too long even on grids 
commonly used by GCMs and ice fiow models today. Most 
of the “intersections” of At and Gj will result in nothing, be- 
cause Ai and Gj are far from each other and obviously do 
not intersect. 

This problem is solved by using an R-tree (Guttman, 1984) 
to avoid having to consider intersections of grid cells that are 
far away from each other. The procedure works as follows: 
first load all the grid cell outlines of A into the R-tree, in- 
dexed by their bounding rectangles. Then loop through each 
grid cell in G, checking the R-tree for any grid cells in A that 
it might intersect with. The polygon intersection algorithm 
is run on each of those grid cell pairs to determine the exact 
outlines of the exchange grid cells. Running time is cut down 
to a more reasonable 0(n logn). 

Past algorithms exist to compute regridding matrices by 
integrating functions around each cell in an exchange grid 
(Ramshaw, 1985; Dukowicz and Kodis, 1987; Jones, 1999). 
These algorithms were originally presented in terms of 
quadrilateral meshes, but they can also be applied to arbi- 
trary meshes. However, they never explicitly compute the ex- 
change grid polygonal outlines. By explicitly computing an 
exchange grid and then using that to produce the regridding 
matrix, we have presented here a procedure that is conceptu- 
ally simpler, possibly more fiexible, but almost certainly not 
faster to run. 


www.geosci-model-dev.net/7/883/2014/ 


Geosci. Model Dev., 7, 883-907, 2014 


906 


R. Fischer et al.: Conservative regridding for ice-atmosphere coupling 


Acknowledgements. Thanks to Ed Bueler, Jeremy Fyeke, William 
Lipseomb, Christian Rodehaeke and Ryan Walker for reading 
drafts. Thanks also to Roger Braithwaite, Nieholas Gould, Regine 
Hoek and David Rind. This researeh is sponsored by the NASA 
MAP program. 

Edited by: D. Roehe 


References 

Asehwanden, A., Bueler, E., Khroulev, C., and Blatter, H.: An en- 
thalpy formulation for glaeiers and iee sheets, J. Glaeiol., 58, 
441^57, doi:10.3189/2012JoGllJ088, 2012. 

Balaji, V., Anderson, J., Held, I., Winton, M., Duraehta, J., Maly- 
shev, S., and Stouffer, R. J.: The Exehange Grid: a meeha- 
nism for data exehange between Earth System eomponents on 
independent grids, in: Parallel Computational Fluid Dynamies: 
Theory and Applieations, Proeeedings of the 2005 International 
Conferenee on Parallel Computational Fluid Dynamies, College 
Park, MD, USA, 24-27 May, edited by: Deane, A., Brenner, G., 
Eeer, A., Emerson, D., MeDonough, J., Periaux, J., Satofuka, N., 
and Tromeur-Dervout, D., Elsevier, 2006. 

Best, M. J., Baljaars, A., Poleher, J., and Viterbo, R: A proposed 
strueture for eoupling tiled surfaees with the planetary boundary 
layer, J. Hydrometerol., 5, 1271-1278, doi:10.1 175/JHM-382.1, 
2004. 

Bindsehadler, R., Nowieki, S., Abe-Ouehi, A., Asehwanden, A., 
Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., 
Herzfeld, U., Jaekson, C., Johnson, J., Khroulev, C., Lev- 
ermann. A., Lipseomb, W, Martin, M., Morlighem, M., 
Parizek, B., Pollard, D., Priee, S., Ren, D., Saito, F., Sato, T., Sed- 
dik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W: 
lee-sheet model sensitivities to environmental foreing and their 
use in projeeting future sea level (the SeaRISE projeet), J. 
Glaeiol., 59, 195-224, doi:10.3189/2013JoG12J125, 2013. 

Box, J. E., Cappelen, J., Chen, C., Deeker, D., Fettweis, X., 
Mote, T., Tedeseo, M., van de Wal, R. S. W, and Wahr, J.: Green- 
land lee Sheet [in Aretie Report Card 2012], available at: http:// 
www.aretie.noaa.gov/reporteard/greenland_iee_sheet.html (last 
aeeess: 3 Deeember 2013), 2013. 

Braden, B.: The surveyor’s area formula. Coll. Math. J., 17, 326- 
337, doi: 10.2307/2686282, 1986. 

Braithwaite, R.: On glaeier energy balanee, ablation, and air tem- 
perature, J. Glaeiol., 27, 381-391, 1981. 

Bueler, E. and Brown, J.: Shallow shelf approximation as a sliding 
law in a thermomeehanieally eoupled iee sheet model, J. Geo- 
phys. Res., 114, F03008, doi: 10.1 029/2008 JFOOl 179, 2009. 

Cgal: Computational Geometry Algorithms Library, available at: 
http://www.egal.org, last aeeess: 3 Deeember 2013. 

Chin, F. and Wang, C. A.: Optimal algorithms for the interseetion 
and the minimum distanee problems between planar polygons, 
IEEE T. Comput., 32, 1203, doi: 1 0. 1 109/TC. 1983. 1676 186, 
1983. 

Chureh, J., Clark, R, Cazenave, A., Gregory, J, Jevrejeva, S., Lev- 
ermann. A., Merrifield, M., Milne, G, Nerem, R., Nunn, R, 
Payne, A., Pfeffer, W, Stammer, D., and Unnikrishnan, A. :Sea 
Level Change. In Climate Change 2013: The Physieal Seienees 
Basis, Contribution of Working Group I to the Fifth Assessment 


Report of the Intergovernmental Panel on Climate Change, Cam- 
bridge University Press, 2013. 

Collins, W. J., Bellouin, N., Doutriaux-Boueher, M., Gedney, N., 
Halloran, R, Hinton, T., Hughes, J., Jones, C. D., Joshi, M., 
Liddieoat, S., Martin, G., O’Connor, R, Rae, J., Senior, C., 
Siteh, S., Totterdell, I., Wiltshire, A., and Woodward, S.: Devel- 
opment and evaluation of an Earth- System model - HadGEM2, 
Geosei. Model Dev., 4, 1051-1075, doi: 10.51 94/gmd-4- 1051- 
2011,2011. 

Comford, S. L., Martin, D. R, Graves, D. T., Ranken, D. R, 
Le Broeq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., 
and Lipseomb, W. H.: Adaptive mesh, finite volume model- 
ing of marine iee sheets, J. Comput. Phys., 232, 529-549, 
doi:10.1016/j.jep.2012.08.037, 2013. 

Dukowiez, J. K. and Kodis, J. W: Aeeurate eonservative remap- 
ping (rezoning) for arbitrary Lagrangian-Eulerian eomputations, 
SIAM J. Sei. Stat. Comp., 8, 305-321, doi: 10.1 137/0908037, 
1987. 

Fiseher, R.: Glint2, available at: http://eitibob.github.io/glint2 (last 
aeeess: 3 Deeember 2013), 2013. 

Gamma, E., Helm, R., Johnson, R., and Vlissides, J.: Design Pat- 
terns: Elements of Reusable Objeet-Oriented Software, 1st Edn., 
Addison- Wesley Professional, 1994. 

Goelzer, H., Huybreehts, R, Furst, J. J., Niek, F. M., Ander- 
sen, M. L., Edwards, T. L., Fettweis, X., Payne, A., J., 
and Shannon, S.: Sensitivity of Greenland iee sheet pro- 
jeetions to model formulations, J. Glaeiol., 59, 733-749, 
doi:10.3189/2013JoG12J182, 2013. 

Gould, N., Orban, D., and Toint, R: GALAHAD, a library 
of thread-safe Fortran 90 paekages for large-seale nonlin- 
ear optimization, ACM T. Math. Software, 29, 353-372, 
dohlO.l 145/962437.962438, 2003. 

Greve, R.: On the response of the Greenland lee Sheet to 
greenhouse elimate ehange, Climatie Change, 46, 289-303, 
doi: 10. 1 023/A: 1 005647226590, 2000. 

Guttman, A.: R- trees - a dynamie index strueture for spatial seareh- 
ing, in: SIGMOD ’84, Proeeedings of the 1984 ACM SIGMOD 
international eonferenee on management of data, 14, 47-57, 
dohlO.l 145/971697.602266, 1984. 

Heath, T. and Arehimedes: The Works of Arehimedes, C. J. Clay 
and Sons, Cambridge University Press Warehouse, available 
at: https://arehive.org/details/worksofarehimede029517mbp (last 
aeeess: 3 Deeember 2013), 1897. 

Hurrell, J. W, Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., 
Kushner, P. J., Lamarque, J. R, Large, W. G., Lawrenee, D., 
Lindsay, K., Lipseomb, W. H., Long, M. C., Mahowald, N., 
Marsh, D. R., Neale, R. B., Raseh, R, Vavrus, S., Vertenstein, M., 
Bader, D., Collins, W. D., Haek, J. J., Kiehl, J., and Mar- 
shall, S.: The Community Earth System Model: a framework for 
eollaborative researeh, B. Am. Meteorol. Soe., 94, 1339-1360, 
doi: 1 0. 1 1 75/BAMS-D- 1 2-00 121.1,2013. 

Huybreehts, R: The present evolution of the Greenland iee sheet: 
an assessment by modelling. Global Planet. Change, 9, 39-51, 
doi: 10.101 6/092 1-8181 (94)90006-X, 1 994. 

Jones, P. W: First- and seeond-order eonservative remap- 
ping sehemes for grids in spherieal eoordinates, Mon. 
Weather Rev., 127, 2204-2210, doi:10.1 175/1520- 

0493(1999)127<2204:FASOCR>2.0.CO;2, 1999. 


Geosei. Model Dev., 7, 883-907, 2014 


www.geosci-model-dev.net/7/883/2014/ 


R. Fischer et al.: Conservative regridding for ice-atmosphere coupling 


907 


Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continen- 
tal scale, high order, high spatial resolution, ice sheet modeling 
using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, 
F01022, doi: 1 0. 1 029/20 11JF002 140, 2012. 

Lauritzen, R H. and Nair, R. D.: Monotone and Conservative 
Cascade Remapping between Spherical Grids (CaRS): Regu- 
lar Latitude-Longitude and Cubed- Sphere Grids, Mon. Weather 
Rev., 136, 1416, dohlO.l 175/2007MWR2181.1, 2008. 

Leung, L. R. and Ghan, S. J.: Parameterizing subgrid oro- 
graphic precipitation and surface cover in climate mod- 
els, Mon. Weather Rev., 126, 3271-3291, doiilO.l 175/1520- 
0493(1998)126<3271:PSOPAS>2.0.CO;2, 1998. 

Lipscomb, W. H., Fyke, J. G., Vizcaino, M., Sacks, W. J., Wolfe, J., 
Vertenstein, M., Craig, A., Kluzek, E., and Lawrence, D. M.: Im- 
plementation and initial evaluation of the Glimmer Community 
Ice Sheet Model in the Community Earth System Model, J. Cli- 
mate, 26, 7352-7371, doi: 10. 1 175/JCLI-D- 12-00557.1, 2013. 

Nowicki, S., Bindschadler, R., Abe-Ouchi, A., Aschwanden, A., 
Bueler, E., Choi, H., Fastook, J., Granzow, G., Greve, R., 
Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Kroulev, C., 
Larour, E., Levermann, A., Lipsomb, W, Martin, M., 

Morlighem, M., Parizek, B., Pollard, D., Price, S., Rignot, E., 
Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Taka- 
hashi, K., Walker, R., and Wang, W: Insights into spatial sen- 
sitivities of ice mass response to environmental change from the 
SeaRISE ice sheet modeling project I: Antarctica, J. Geophys. 
Res., 118, 1002-1024, doi: 10. 1002/jgrf.20081, 2013a. 

Nowicki, S., Bindschadler, R., Abe-Ouchi, A., Aschwanden, A., 
Bueler, E., Choi, H., Fastook, J., Granzow, G., Greve, R., 
Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Kroulev, C., 
Larour, E., Levermann, A., Lipsomb, W, Martin, M., 

Morlighem, M., Parizek, B., Pollard, D., Price, S., Rignot, E., 
Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Taka- 
hashi, K., Walker, R., and Wang, W: Insights into spatial sen- 
sitivities of ice mass response to environmental change from the 
SeaRISE ice sheet modeling project II: Greenland, J. Geophys. 
Res.-Earth, 118, 1025-1044, doi:10.1002/jgrf.20076, 2013b. 

Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheet- 
shelf model, and application to Antarctica, Geosci. Model Dev., 
5, 1273-1295, doi: 1 0.5 194/gmd-5- 1273-20 12, 2012. 

Qu, X. and Hall, A.: Assessing snow albedo feedback in 
simulated climate change, J. Climate, 19, 2617-2630, 
doi:10.1175/JCLI3750.1, 2006. 

Ramshaw, J. D.: Conservative rezoning algorithm for general- 
ized two-dimensional meshes, J. Comput. Phys., 59, 193-199, 
doi: 10. 101 6/002 1-9991 (85)90141 -X, 1985. 

Ridley, J. K., Huybrechts, P, Gregory, J. M., and Lowe, J. A.: Elim- 
ination of the Greenland Ice Sheet in a high CO 2 climate, J. Cli- 
mate, 18, 3409-3427, doi:10.1175/JCLI3482.1, 2005. 


Ridley, J., Gregory, J. M., Huybrechts, P, and Lowe, J.: Thresholds 
for irreversible decline of the Greenland ice sheet, Clim. Dynam., 
35, 1059-1057, doi:10.1007/s00382-009-0646-0, 2010. 

Rutt, I. C., Hagdom, M., Hulton, N. R. J., and Payne, A. J.: The 
Glimmer community ice sheet model, J. Geophys. Res., 114, 
F02004, doi:10.1029/2008JF001015, 2009. 

Schmidt, G. A., Bitz, C. M., Mikolajewicz, U., and Tremblay, L. B.: 
Ice-ocean boundary conditions for coupled models. Ocean 
Model., 7, 59-74, doi:10.1016/S1463-5003(03)00030-l, 2004. 

Schmidt, G. A., Ruedy, R., Hansen, J. E., Aleinov, I., Bell, N., 
Bauer, M., Bauer, S., Cairns, B., Canute, V, Cheng, Y., Del Ge- 
nio. A., Faluvegi, G., Friend, A. D., Hall, T. M., Hu, Y, Kel- 
ley, M., Kiang, N. Y, Koch, D., Lacis, A. A., Lemer, J., Lo, K. K., 
Miller, R. L., Nazarenko, L., Oinas, V, Perlwitz, J. P, Perl- 
witz, J., Rind, D., Romanou, A., Russell, G. L., Sato, M., Shin- 
dell, D. T., Stone, P. H., Sun, S., Tausnev, N., Thresher, D., and 
Yao, M.-S.: Present day atmospheric simulations using GISS 
ModelE: comparison to in-situ, satellite and reanalysis data, J. 
Climate, 19, 153-192, doi: 1 0. 1 175/JCLI36 12.1, 2006. 

Snyder, J. P: Map projections - a working manual, available at: 
http://pubs.er.usgs.gov/publication/ppl395 (last access: 3 De- 
cember 2013), 1987. 

Stone, E. J., Lunt, D. J., Rutt, I. C., and Hanna, E.: Investigating the 
sensitivity of numerical model simulations of the modern state of 
the Greenland ice-sheet and its future response to climate change. 
The Cryosphere, 4, 397-417, doi:10.5194/tc-4-397-2010, 2010. 

van den Broeke, M., Smeets, P, Ettema, J., and Munneke, P. K.: 
Surface radiation balance in the ablation zone of the 
west Greenland ice sheet, J. Geophys. Res., 113, D13105, 
doi:10.1029/2007JD009283, 2008. 

Vaughan, D., Comiso, J., Allison, I., Carrasco, J., Kaser, G., Kwok, 
R., Mote, P, Murray, T., Paul, F., Ren, J., Rignot, E., Solom- 
ina, O., Steffen, K. and Zhang, T.: Observations: Cryosphere. In 
Climate Change 2013: The Physical Sciences Basis, Contribu- 
tion of Working Group I to the Fifth Assessment Report of the 
Intergovernmental Panel on Climate Change, Cambridge Univer- 
sity Press, 2013. 

Vizcaino, M., Lipscomb, W, Sacks, W, van Angelen, J. H., 
Wouters, B., and van den Broeke, M. R.: Greenland surface mass 
balance as simulated by the Community Earth System Model. 
Part I: model validation and 1850-2005 results, J. Climate, 26, 
7793-78 1 2, doi: 1 0. 1 1 75/JCLI-D- 1 2-006 15.1,2013. 

Zagorodnov, V, Nagornov, O., and Thompson, L. G.: Influ- 
ence of air temperature on a glacier’s active-layer temperature, 
Ann. Glaciol., 43, 285-291, doi: 10.3 189/172756406781812203, 
2006. 

Zienkiewicz, O. C., Taylor, R. L., and Nithiarasu, R: The Finite 
Element Method for Fluid Dynamics, 7th Edn., Butterworth- 
Heinemann, 2013. 


www.geosci-model-dev.net/7/883/2014/ 


Geosci. Model Dev., 7, 883-907, 2014 


